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ABSTRACT 

We study a solid protoplanetary core undergoing radial migration in a protoplanetary disk. We 
consider cores in the mass range ~ 1 — 10 M ffi embedded in a gaseous protoplanetary disk at different 
radial locations. 

We suppose the core luminosity is generated as a result of planetesimal accretion and calculate the 
structure of the gaseous envelope assuming hydrostatic and thermal equilibrium. This is a good approx- 
imation during the early growth of the core while its mass is less than the critical value, M cr i t , above 
which such static solutions can no longer be obtained and rapid gas accretion begins. The critical value 
corresponds to the crossover mass above which rapid gas accretion begins in time dependent calculations. 

We model the structure and evolution of the protoplanetary nebula as an accretion disk with constant 
a. We present analytic fits for the steady state relation between disk surface density and mass accretion 
rate as a function of radius. 

We calculate M cr i t as a function of radial location, gas accretion rate through the disk, and planetesimal 
accretion rate onto the core. For a fixed planetesimal accretion rate, M cr u is found to increase inwards. 
On the other hand it decreases with the planetesimal accretion rate and hence the core luminosity. 

We consider the planetesimal accretion rate onto cores migrating inwards in a characteristic time 
~ 10 3 — 10 5 yr at 1 AU as indicated by recent theoretical calculations. We find that the accretion rate 
is expected to be sufficient to prevent the attainment of M cr i t during the migration process if the core 
starts off significantly below it. Only at those small radii where local conditions are such that dust, and 
accordingly planetesimals, no longer exist can M cr i t be attained. 

At small radii, the runaway gas accretion phase may become longer than the disk lifetime if the mass 
of the core is too small. However, within the context of our disk models, and if it is supposed that some 
process halts the migration, massive cores can be built-up through the merger of additional incoming 
cores on a timescale shorter than for in situ formation. A rapid gas accretion phase may thus begin 
without an earlier prolonged phase in which planetesimal accretion occurs at a reduced rate because of 
feeding zone depletion in the neighborhood of a fixed orbit. 

Accordingly, we suggest that giant planets may begin to form through the above processes early in 
the life of the protostellar disk at small radii, on a timescale that may be significantly shorter than that 
derived for in situ formation. 

Subject headings: accretion, accretion disks — solar system: formation — planetary systems 



gravitational potential energy liberated by matter flow- 
ing inwards at a rate M ~ 10 _8±1 Mq yr -1 (Hartmann 
et al. 1998). The nonobservation of disks around older 
T Tauri stars together with these values of M suggest a 
disk lifetime of between 10 6 and 10 7 yr (Strom, Edwards & 
Skrutskie 1993). Masses between 10~ 3 and 10 _1 M and 
dimensions in the range 10-100 AU have been estimated 
(Beckwith & Sargent 1996). 

Most theoretical protostellar disk models have relied on 
the a-parametrization proposed by Shakura & Sunyaev 
(1973). In this context, the disk anomalous turbulent 
viscosity, which enables angular momentum to be trans- 
ported outwards and therefore matter to flow inwards and 
be ultimately accreted by the central star, is assumed to 
give rise to a stress tensor which is simply proportional to 
the gas pressure. So far, only MHD instabilities (Balbus & 
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1. INTRODUCTION 

The hypothesis that the planets in the solar system 
were formed in a flattened differentially rotating gaseous 
disk was originally proposed by Kant (1755) and Laplace 
(1796). Since then, the presence of disks around low-mass 
young stellar objects has been inferred from their infrared 
excess (Adams, Lada & Shu 1987). Recently they have 
also been imaged directly with the Hubble Space Telescope 
(McCaughrean & O'Dell 1996; Burrows et al. 1996; Mc- 
Caughrean et al. 1998; Krist et al. 1998; Stapelfeldt et 
al. 1998). Surveys in the Orion nebula (Stauffer et al. 
1994) and the Taurus- Auriga dark clouds (Beckwith et al. 
1990) indicate that these disks are common, apparently 
surrounding between 25 and 75% of the young stellar ob- 
jects. Their infrared emission may be produced by the 
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Hawley 1991) have been shown to be able to produce and 
sustain turbulence in accretion disks, and they do lead 
to an a-type disk (Balbus & Hawley 1998). However, be- 
cause these instabilities develop only in an adequately ion- 
ized fluid, they may not operate everywhere in protostellar 
disks (Gammie 1996). Therefore it is likely that the pa- 
rameter a, to which the viscosity is simply related, is not 
constant through these disks. It may even be that only 
parts of these disks can be described using this a prescrip- 
tion for the viscosity. However, we may still learn about 
disks from these models in the same way as we learned 
about stars from simple polytropic models. Therefore, for 
the purpose of considering planet formation, as we are in- 
terested in here, we will use such models. 

Planets are believed to form out of protostellar disks 
by either gravitational instability (Kuiper 1951; Cameron 
1978; Boss 1998) or by a process of growth through plan- 
etesimal accumulation followed, in the giant planet case, 
by gas accretion (Safronov 1969; Wetherill & Stewart 1989; 
Perri & Cameron 1974; Mizuno 1980; Bodenheimer & Pol- 
lack 1986). The first mechanism is expected to produce 
preferentially massive objects in the outer parts of the 
disk, if anything. In this paper we will study planetary 
formation within the context of the second mechanism, 
which is commonly accepted as the most likely process by 
which planets form in at least the inner ten astronomical 
units of protostellar disks. We note however that impor- 
tant issues related to this model still remain to be resolved 
(see Lissauer 1993 for a review). 

Up to now, planetary formation has been studied at 
a given location in a disk. Most work has concentrated 
on orbital distances corresponding to the neighborhood 
of Jupiter. More recently, prompted by the detection of 
planets orbiting at short distances from their host star, in 
situ formation of giant planets at these locations has also 
been considered (Ward 1997a; Bodenheimer 1998; Boden- 
heimer, Hubickyj & Lissauer 1998). 

However, the importance of orbital migration as recently 
indicated by both observations (see Marcy and Butler 1998 
and references therein) and theory (see Lin & Papaloizou 
1993 and references therein; Ward 1997b) suggests that 
planets may not form at a fixed location in the disk, but 
more likely grow while migrating through the nebula. It 
is the purpose of this paper to investigate the effect of 
migration on planetary formation. 

In § H we construct steady state a-disk models for 
a = 10 and a = 10 -3 . A range of gas accretion rates, 
M, varying between 10~ 6 and 10~ 9 M Q yr -1 are consid- 
ered. The steady state assumption is reasonable in the 
inner regions, below 5 AU from the central star, where the 
local viscous timescale is short, typically on the order of 
10 4 yr for M = 10~ 7 M Q yr -1 . We present analytic fits for 
the steady state relation between disk surface density and 
mass accretion rate as a function of disk location. These 
fits can be used to solve the diffusion equation which gov- 
erns the disk evolution. 

In § U we describe the construction of the protoplanet 
models based on a solid core with gaseous enve lope . We re- 
view the theory of giant planet formation in § 3.1 and give 
the e qua tions governing a protoplanet atmosphere in § |3.2 
In 



|3.3| we give the results of numerical calculations for 
the critical core mass above which the atmosphere cannot 



remain in hydrostatic and thermal equilibrium but must 
evolve, with the protoplanet entering a rapid gas accretion 
phase. These calculations are done at different locations in 
our disk models. For a given planetesimal accretion rate 
which supplies the core luminosity necessary to support 
the envelope, critical core masses are found to increase sig- 
nificantly by factors of 2-3 between 5 AU and 0.05 AU, at 
which point local conditions may not enable planetesimals 
to exist. 

In § |] we consider migration of the protoplanetary cores 
which, according to recent estimates of the effects of tidal 
interactions with the disk by Ward (1997b), may occur 
on a timescale between 10 4 yr and 10 5 yr for a core mass 
of several earth masses at 5 AU. This is also comparable 
to the proposed formation time. We perform simulations 
which indicate that such a migrating protoplanet is likely 
to accrete 24% or more of the planetesimals initially in- 
terior to it. This accretion is likely to maintain the core 
luminosity such that attainment of the critical mass does 
not occur until it reaches small radii ~ 0.1 AU where plan- 
etesimals no longer exist. 
In 
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we discuss our results in the context of short 
period giant planets. We point out that the processes in- 
vestigated in this paper are likely to result in a giant planet 
orbiting at small radii on a timescale significantly shorter 
than that derived for in situ formation. Finally in § 5.2 
we summarize our results. 

2. DISK MODELS 

2.1. Vertical Structure 

2.1.1. Basic Equations 

Here we consider the equations governing the disk verti- 
cal structure in the thin-disk approximation. Using cylin- 
drical coordinates (r, ip, z) based on the central star and 
such that the z = plane corresponds to the disk mid- 
plane, we adopt the equation of vertical hydrostatic equi- 
librium in the form: 



1 dP 

p dz 



-n 2 z, 



(i) 



together with the energy equation, which states that the 
rate of energy removal by radiation is locally balanced by 
the rate of energy production by viscous dissipation: 
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(2) 



where F is the radiative flux of energy through a surface 
of constant z which is given by: 



F 



-16oT 3 dT 

3np dz 



(3) 



Here p is the mass density per unit volume, P is the 
pressure, T is the temperature, fl is the angular veloc- 
ity, v is the kinematic viscosity, k is the opacity, which in 
general depends on both p and T, and a is the Stefan- 
Boltzmann constant. We consider thin disks which are 
in Keplerian rotation around a star of mass M», so that 
Q 2 = GA'F/r 3 , G being the gravitational constant. In 
writing equation (||), we have assumed that the disk at 
the radius considered is optically thick. However, when 
the disk is optically thin, i.e. when up integrated over the 
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disk thickness is small compared to unity, the temperature 
gradient given by equation (|J) is small, so that the results 
we get are consistent in that case also. 

To close the system of equations, we relate P, p and T 
through the equation of state of an ideal gas: 



P = 



pkT 
pm H : 



(4) 



where k is the Boltzmann constant, p, is the mean molec- 
ular weight and tuh is the mass of the hydrogen atom. 
Here we shall limit our calculations to temperatures lower 
than 4,000 K, so that, at the densities of interest, hydro- 
gen is in molecular form. Since the main component of 
protostellar disks is hydrogen, it is a reasonable approxi- 
mation to take p — 2. Tests using a more sophisticated 
equation of state such as that of Chabrier et. al. (1992) 
indicate only minor differences in the results for the range 
of temperatures and densities considered. Similarly trans- 
port of energy by convection can be neglected here (see 
also Lin & Papaloizou 1985). 

We adopt the a-parametrization of Shakura & Sun- 
yaev (1973), so that the kinematic viscosity is written 
v = ac^/ri, where c s is the isothermal sound speed 
(Cg = P/p). Although in general a may be a function 
of both r and z, we shall limit our calculations presented 
below to cases with constant a (see discussion in § |J). 
With this formalism, equation (0) becomes: 



dF 9 n 
2.1.2. Boundary Conditions 



(5) 



We have to solve three first order ordinary differential 
equations for the three variables F, P (or equivalently p), 
and T as a function of z at a given radius r. Accordingly, 
we need three boundary conditions at each r. We denote 
with a subscript s values at the disk surface. 

A boundary condition is obtained by integrating equa- 
tion (||) over z between —H and H which define the lower 
and upper boundary of the disk, respectively. Since by 
symmetry F(z = 0) = 0, this gives: 



-MM 2 , 



(6) 



where we have defined M st = 37r(z^)£ with £ = J H H pdz 

being the disk surface mass density and (is) = J^ H pvdz/T, 
being the vertically averaged viscosity. If the disk were in 
a steady state, M st would not vary with r and would be 
the constant accretion rate through the disk. In ge nera l 
however, the disk is not in a steady state (but see § |2.2| ) 
such that this quantity does depend on r and the disk 
undergoes time dependent evolution. 

Another boundary condition is obtained by integrating 
equation (jl]) over z between H and infinity. A detailed 
derivation of this condition is presented in Appendix ^j. 
Here we just give the result: 



Ps 



(7) 



where r a b is the optical depth above the disk. This con- 
dition is familiar in stellar structure, where fl 2 H would 



be replaced by the acceleration of gravity at the stellar 
surface (e.g. Schwarzschild 1958). Since we have defined 
the disk surface such that the atmosphere above the disk 
is isothermal, we have to take T a b <C 1. Providing this is 
satisfied, the result s do not depend on the value of T a b we 
choose (see § |2.1.3| ). 

A third and final boundary condition is given by the ex- 
pression of the surface temperature (see Appendix [A| for a 
detailed derivation of this expression): 



2a IT: -Tu 
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—M, f n z 



0. 



(8) 
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Here the disk is assumed immersed in a medium with 
background temperature T},. The surface opacity k s in 
general depends on both T s and p s and we have used 
c 2 s = kT '/ \fimjj). The boundary condition (|^) is the same 
as that used by Levermore & Pomraning (1981) in the Ed- 
dington approximation (their eq. [56] with 7=1/2). In 
the simple case when Tf, — and the surface dissipation 
term involving a is set to zero, with M st being retained, it 
simply relates the disk surface temperature to the emer- 
gent radiation flux. 

2.1.3. Model Calculations 

At a given radius r and for a given value of the param- 
eters M st and a, we solve equations (|TJ) (||) and (||) with 
the boundary conditions (^), (Q) and (g) to find the de- 
pendence of the state variables on z. The opacity is taken 
from Bell & Lin (1994). This has contributions from dust 
grains, molecules, atoms and ions. It is written in the form 

p a T b where Kj, a and b vary with temperature. 
The equations are integrated using a fifth-order Runge- 
Kutta method with adaptive step length (Press et al. 
1992). For a specified M st and a, we first calculate the 
surface flux F s and temperature T s from equations (^) 
and (||), respectively. If T s is smaller than about 1,000 K, 
the opacity k s at the disk surface does not depend on p s . 
For larger values of T s , it turns out that the second term 
in equation (||), which contains k s , is negligible compared 
to the third term. Therefore T s can always be calculated 
independently of p s (or equivalently P s ). We then deter- 
mine the value of H, the vertical height of the disk surface, 
iteratively. Starting from an estimated value of H, we cal- 
culate P s and integrate the equations from H down to the 
midplanc z = 0. The condition that F = at z = will 
not in general be satisfied. An iteration procedure is then 
used to adjust the value of H until F = 0atz=0toa 
specified accuracy. 

An important point to note is that as well as finding the 
disk structure, we also determine the surface density £ for 
a given M st = 3tt(z/)£. In this way, a relation between (v) 
and £ is derived. 

In the calculations presented here, we have taken M* = 
1 Mq, the optical depth of the atmosphere above the 
disk surface r a (, = 10~ 2 and a background temperature 
Tb = 10 K. In the optically thick regions of the disk, the 
value of H is independent of the value of T a b we choose. 
However, this is not the case in optically thin regions where 
we find that, as expected, the smaller r a b, the larger H. 
However, this dependence of H on r a b has no physical 
significance, since the surface mass density, the optical 
thickness through the disk and the midplane temperature 
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hardly vary with r a b- This is because the mass is concen- 
trated towards the disk midplane in a layer with thickness 
independent of T a b- 

2.2. Time Dependent Evolution and Quasi-Steady States 

In general an accretion disk is not in a steady state but 
undergoes time dependent evolution. The global evolution 
of the disk is governed by the well known diffusion equa- 
tion for the surface density which can be written in the 
form (see Lynden-Bell & Pringle 1974, Papaloizou & Lin 
1995 and references therein): 



as 

~dt 



3d_ 

r dr 



-1/2 



d_ 

Or 



/2 



(9) 



The characteristic diffusion time scale at radius r is then: 



i{v) 3a 



i_(-V 

la \HJ 



(10) 



For disks with approximately constant aspect ratio H/r, as 
applies to the models considered here, i„ scales as the local 
orbital period. One thus expects that the inner regions re- 
lax relatively quickly to a quasi-steady state which adjusts 
its accretion rate according to the more slowly evolving 
outer parts (see Lynden-Bell & Pringle 1974 and Lin & 
Papaloizou 1985). For estimated sizes of protostellar disks 
of about 50 AU (Beckwith & Sargent 1996), the evolution- 
ary timescale associated with the outer parts is about 30 
times longer than that associated with the inner parts with 
r < 5 AU, which we consider here in the context of plan- 
etary formation. Thus these inner regions are expected to 
be in a quasi-steady state through most of the disk life- 
time. We have verified that this is the case by considering 
solutions of equation (|9|) . 

In order to investigate these solutions and for other pur- 
poses we found it convenient to utilize analytic piece-wise 
p ower law fits to the relation derived above (see 

2.1.3). Details of these fits are given in Appendix |b| 



In Figures |l|a— b we plot both the curves M st (£) obtained 
from the vertical structure integrations and those obtained 
from the piece-wise power law fits. Figures |I|a and [l]b are 
for a = 10 -2 and a = 10~ 3 , respectively. In each case 
the radius varies between 0.01 and 100 AU, and the cal- 
culatio ns are limited to temperatures lower than 4,000 K 
(see § |2.1.1 ). The average errors are 18% and 13% for 
a = 10~ J and 10~ 2 , respectively. Thus the fits give an 
adequate approximation. 

2.3. Steady State Models 

Here, we present solutions corresponding to steady state 
accretion disks. In Figures ||a-c and ||a-c, we plot H/r, 
E and the midplane temperature T m versus r for M st be- 
tween 10 -9 and 10~ 6 M Q yr _1 (assuming this quantity is 
the same at all radii, i.e. the disk is in a steady state) for 
a = 10~ 2 and a — 10~ 3 , respectively. 

An inspection of Figures |a and ||a indicates that the 
outer parts of the disk are shielded from the radiation of 
the central star by the inner parts, apart possibly from 
the very outer parts, which are optically thin anyway and 
therefore do not reprocess any radiation. This is in agree- 
ment with the results of Lin & Papaloizou (1980) and 
Bell et al. (1997). For a = 10~ 3 , the radius beyond 



which the disk is not illuminated by the central star varies 
from 0.2 AU to about 3 AU when M st goes from 10~ 9 
to 10 -6 Mq yr -1 . These values of the radius move to 
0.1 and 2 AU when a = 10~ 2 . Since reprocessing of the 
stellar radiation by the disk is not an important heating 
factor below these radii, this process will in general not be 
important in these models of protostellar disks. We note 
that this result is independent of the value of t q {, we have 
taken. Indeed, as we pointed out above, only the thickness 
of the optically thin parts of the disk gets larger when T a & 
is decreased. 

However, there are some indications that disks such as 
HH30 may be flared (Burrows et al. 1996). In this context, 
Chiang & Goldreich (1997) have considered a model based 
on reprocessing in which the dust and gas are at different 
temperatures. However, as they pointed out, some issues 
regarding this model remain to be resolved. In any case, 
it is possible that a multiplicity of solutions exists when 
reprocessing is taken into account, with it being important 
for cases in which the disk is flared and unimportant when 
it is not, such as maybe HK Tau (Stapelfeldt et al. 1998, 
Koresko 1998). 

The values of H/r, £ and T m we get are similar to those 
obtained by Lin & Papaloizou (1980), who adopted a pre- 
scription for viscosity based on convection, and Bell et al. 
(1997). Since H is measured from the disk midplane to 
the surface such that T a b is small, it is larger than what 
would be obtained if a value of 2/3 were adopted for r a b, 
as is usually the case. However, this does not affect other 
physical quantities. We also recall that H, as defined here, 
is about 2-3 times larger than c s /f2, with c s being the mid- 
plane sound speed, which is commonly used to define the 
disk semithickness. 

In this paper we shall consider the migration of proto- 
planetary cores from ~ 1-5 AU, where they are supposed 
to form under conditions where ice exists, down to the 
disk inner radii (see § 0) . It is therefore of interest to esti- 
mate the mass of planetesimals contained inside the orbit 
of a core when it forms, since this can potentially be ac- 
creted by the core during its migration. From Figures ||b 
and we estimate the mass of planetesimals M p (r) con- 
tained within a radius r using M p (r) = 10~ 2 x 7rr 2 E(r) 
for r = 1 AU and r = 5 AU. Here we have assumed 
a gas to dust ratio of 100. The values of M p corre- 
sponding to a = 10~ 2 and 10~ 3 and M = 10~ 6 and 
10~ 7 Mq yr" 1 are listed in Table [l[ We have checked 

that the disk with a = 10~ 3 and M = 10~ 6 Mq yr _1 , al- 
though relatively massive, is gravitationally stable locally. 
Namely the Toomre parameter Q = (M*/M(r))(H/r), 
with M(r) — 7rr 2 S(r), is larger than unity. 

It is also of interest, in relation to the possibility of gi- 
ant planets being located at small radii, to estimate the 
mass of gas contained within a radius of 0.1 AU. Figure 0b 
indicates that, when a = 10~ 3 , this mass is about 073 
Jupiter mass for M > 10~ 7 M Q yr _1 . For a — 10~ 2 , 
there is a similar mass of gas for M > 10~ 6 M Q yr _1 
(see Figure ||d). For typical mass throughput of about 
10 _2 -10 _1 M Q , the lifetime of such a state can range be- 
tween 10 4 and 10 6 yr. Supposing the disk to be termi- 
nated at some small inner radius, this suggests that, if a 
suitable core can migrate there, it could accrete enough 
gas to become a giant planet within the disk lifetime. We 
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note however that the conditions for that to happen are 
marginal even in the early stages of the life of the disk 
when M > 1CT 6 - 1CT 7 M© yr" 1 . 

At later stages, when M ~ 10~ 8 M Q yr -1 , the mod- 
els resemble conditions expected to apply to the mini- 
mum mass solar nebula with S ~ 200 g cm~ 2 at 5 AU 
if a = 10~ 2 . Under these conditions, the mass of gas at 
r < 0.1 AU is between 1 and 9 M® for a between 10 -2 
and 10~ 3 . 

Lin et al. (1996) suppose that the inner disk is termi- 
nated by a magnetospheric cavity. In this case migration 
might be supposed to cease if the core is sufficiently far in- 
side it. But gas accretion is likely to be very much reduced 
in this case. 

However, the protoplanet may be able to accrete more 
gas than the mass estimated above if circumstances were 
such that the core migration is halted at some small ra- 
dius before the disk is terminated. We note that Ward 
(1986) finds the direction of type I migration is insensitive 
to the disk surface density profile but that it could reverse 
from inwards to outwards if the disk midplane tempera- 
ture decreased inwards faster than approximately linearly. 
Such a condition would not be expected in the disk models 
considered here. 

On the other hand, conditions may be very different 
if interaction with a stellar magnetic field becomes im- 
portant. It is expected that this happens when magnetic 
and viscous torques become comparable (e.g. Ghosh & 
Lamb 1979). In this region there may be open field lines 
connected to the disk with an outflowing wind (Paatz & 
Camenzind 1996). Such a wind may provide an additional 
angular momentum and energy loss mechanism for the disk 
material (Papaloizou & Lin 1995). The inner regions could 
then be cooler than expected from the constant a models 
considered here leading to a reversal of type I migration. A 
faster gas inflow rate may also prevent the onset of type II 
migration. Thus, although details are unclear, continued 
accretion of inflowing disk gas by a protoplanetary core 
that has stopped migrating in the inner disk may be pos- 
sible. 

3. PROTOPLANETARY CORE GROWTH AND 
EQUILIBRIUM ENVELOPE 

3.1. Background: Formation of Giant Planets 

The solid cores of giant planets are believed to be formed 
via solid body accretion of km-sized planetesimals, which 
themselves are produced as a result of the sedimentation 
and collisional growth of dust grains in the protoplanetary 
disk (see Lissauer 1993 and references therein). Once the 
solid core becomes massive enough to gravitationally bind 
the gas in which it is embedded (typically at a tenth of an 
earth mass), a gaseous envelope begins to form around the 
core. 

The build-up of the atmosphere has first been consid- 
ered in the context of the so-called 'core instability' model 
by Perri & Cameron (1974) and Mizuno (1980; see also 
Stevenson 1982 and Wuchterl 1995). In this model, the 
solid core grows in mass along with the atmosphere in 
quasi-static and thermal equilibrium until the core reaches 
the so-called 'critical core mass' above which no equilib- 
rium solution can be found for the atmosphere. As long 
as the core mass is smaller than the critical core mass, 



the energy radiated from the envelope into the surround- 
ing nebula is compensated for by the gravitational energy 
which the planetesimals entering the atmosphere release 
when they collide with the surface of the core. During this 
phase of the evolution, both the core and the atmosphere 
grow in mass relatively slowly. By the time the core mass 
reaches the critical core mass, the atmosphere has grown 
massive enough so that its energy losses can no longer be 
compensated for by the accretion of planetesimals alone. 
At that point the envelope has to contract gravitationally 
to supply more energy. This is a runaway process, leading 
to the very rapid accretion of gas onto the protoplanet and 
to the formation of giant planets such as Jupiter. In ear- 
lier studies it was assumed that this rapid evolution was a 
dynamical collapse, hence the designation 'core instability' 
for this model. 

Further time-dependent numerical calculations of pro- 
toplanetary evolution by Bodcnhcimer & Pollack (1986) 
support this model, although they show that the core mass 
beyond which runaway gas accretion occurs, which is re- 
ferred to as the 'crossover mass', is slightly larger than the 
critical core mass, and that the very rapid gravitational 
contraction of the envelope is not a dynamical collapse. 
The designation 'crossover mass' comes from the fact that 
rapid contraction of the atmosphere occurs when the mass 
of the atmosphere is comparable to that of the core. Once 
the crossover mass is reached, the core no longer grows 
significantly. 

More recent simulations by Pollack et al. (1996) show 
that the evolution of a protoplanet is governed by three 
distinct phases. During phase 1, runaway planetesimal ac- 
cretion occurs which leads to the depletion of the feeding 
zone of the protoplanet. At this point, when phase 2 be- 
gins, the atmosphere is massive enough that the location 
of its outer boundary is determined by both the mass of 
gas and planetesimals. As more gas is accreted, this outer 
radius moves out, so that the feeding zone is increased and 
more planetesimals can be captured, which in turn enables 
more gas to enter the atmosphere. The protoplanet grows 
in this way until the core reaches the crossover mass, at 
which point runaway gas accretion occurs and phase 3 be- 
gins. The timescale for planet formation is determined al- 
most entirely by phase 2, and is found to be a few million 
years at 5 AU. For typical disk models, this is comparable 
to the disk lifetime. Note however that isolation of the pro- 
toplanetary core, as it occurs at the end of phase 1, may be 
prevented under some circumstances by tidal interaction 
with the surrounding gaseous disk (Ward & Hahn 1995). 

Conditions appropriate to Jupiter's present orbital ra- 
dius are normally considered and then the critical or 
crossover mass is found to be around 15 M®. This is con- 
sistent with models of Jupiter which indicate that it has a 
solid core of about 5-15 M© (Podolak et al. 1993). 

We note that although phase 3 is relatively rapid com- 
pared to phase 2 for conditions appropriate to Jupiter's 
present location, it may become longer when the luminos- 
ity provided by the accretion of planetesimals, and hence 
the cr itica l core mass, is reduced (see Pollack et al. 1996 
and § 5.1). The designation 'runaway' or 'rapid' gas ac- 
cretion may then become confusing. 

The similarity between the critical and crossover masses 
is due to the fact that, although there is some liberation 
of gravitational energy as the atmosphere grows in mass 



6 



together with the core, the effect is small as long as the 
atmospheric mass is small compared to that of the core. 
Consequently the hydrostatic and thermal equilibrium ap- 
proximation for the atmosphere is a good one for core 
masses smaller than the critical value. Therefore we use 
this approximation here and investigate how the critical 
core mass varies with location and physical conditions in 
the protoplanetary disk. 

3.2. Basic Equations Governing a Protoplanetary 
Envelope 

Let R be the spherical polar radius in a frame with ori- 
gin at the center of the protoplanet's core. We assume 
that we can model the protoplanet as a spherically sym- 
metric nonrotating object. We also assume that it is in 
hydrostatic and thermal equilibrium. The equation of hy- 
drostatic equilibrium is then: 



dP 
dR 



-9P, 



(11) 



where g = GM/R 2 is the acceleration due to gravity, 
M(R) being the mass contained in the sphere of radius 
R (this includes the core mass if R is larger than the core 
radius). We also have the definition of density: 



dM 
~dR 



4irR 2 p. 



(12) 



At the high densities that occur at the base of a pro- 
toplanetary envelope, the gas cannot be considered to be 
ideal. Thus we adopt the equation of state for a hydrogen 
and helium mixture given by Chabrier et al. (1992). We 
adopt the mass fractions of hydrogen and helium to be 0.7 
and 0.28, respectively. We also use the standard equation 
of radiative transport in the form: 



dT 
dR 



-3np 



16crT 3 4nR 2 ' 



(13) 



Here L is the radiative luminosity. Denoting the radiative 
and adiabatic temperature gradients by V ra d and V a( j, re- 
spectively, we have: 



<91nT 



1 = 



3hiL core P 



dlnPJ rad QAnaGMT^ 



and 



<91nT 
<91nP 



(14) 



(15) 



with the subscript s denoting evaluation at constant en- 
tropy. 

We assume that the only energy source comes from the 
accretion of planctesimals onto the core which, as a result, 
outputs a total core luminosity L core , given by: 



GM core M t 



core ±v± core 



(16) 



Here M core and r core are respectively the mass and the 
radius of the core, and M core is the planetesimal accre- 
tion rate. The luminosity L core is supplied by the gravita- 
tional energy which the planetesimals entering the planet 



atmosphere release near the surface of the core (see, e.g., 
Mizuno 1980; Bodenheimer & Pollack 1986). 
If V ra( 2 < V d, there is stability to convection and thus all 
the energy is transported by radiation, i.e. L = L core . 
When V ra d > V d, there is instability to convection. 
Then, part of the energy is transported by convection, and 
Lcore = L+L conv , where L conv is the luminosity associated 
with convection. We use mixing length theory to evaluate 
Lconv (Cox & Giuli 1968). Then 



irR 2 C p K 2 ml 



dR) \dRJ 



3/2 



:P9 



dp 

dT 



(17) 



where A m ; = \a m iP / {dP / dR)\ is the mixing length, 
a m i being a constant of order unity, (8T/dR) s = 
V a dT (d\n P / dR) , and the subscript P means that the 
derivative has to be evaluated for a constant pressure. 
All the required thermodynamic parameters are given by 
Chabrier et al. (1992). In the numerical calculations pre- 
sented below we fix a m i = 1. 

3.2.1. Inner Boundary 

We suppose that the planet core has a uniform mass 
density p CO re- The composition of the planetesimals and 
the high temperatures and pressures at the surface of the 
core suggest p CO re = 3.2 g cm~ 3 (see Bodenheimer & Pol- 
lack 1986 and Pollack et al. 1996), which is the value we 
will adopt throughout. The core radius, which is the inner 
boundary of the atmosphere, is then given by: 



3M r , 



1/3 



At R : 



^Pco; 

the total mass is equal to M c , 



(18) 



3.2.2. Outer Boundary 

We take the outer boundary of the atmosphere to be at 
the Roche lobe radius r^ of the protoplanet. Thus: 



rr = 2 (MlL 
3 \3M* 



1/3 



(19) 



where M p i = M core + M atm is the planet mass, M atm be- 
ing the mass of the atmosphere, and r is the orbital radius 
of the protoplanet in the disk. 

To avoid confusion, we will denote the disk midplane 
temperature, pressure and mass density at the distance r 
from the central star by T m , P m and p m , respectively. 



At R 
to P m 



= tl, the mass is equal to M p i, the pressure is equal 
and the temperature is given by 



T = 



4 3t l L 

in. ' 



IQitarj^ 



1/4 



(20) 



where we approximate the additional optical depth above 
the protoplanet atmosphere, through which radiation 
passes, by: 



tl = n(p rn ,T m ) p m r L . 



(21) 
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3.3. Calculations 

For a particular disk model, at a chosen radius r, for 
a given core mass M core and planetesimal accretion rate 
M core , we solve the equations (|TT|), (12) and ( |l5| ) with the 
boundary conditions described above to get the structure 
of the envelope. The opacity law adopted is the same as 
that for the disk models. In general, the deep interior 
of the envelope becomes convective with the consequence 
that the value of the opacity does not matter there. 

For a fixed M core at a given radius, there is a critical 
core mass M cr u above which no solution can be found, 
i.e. there can be no atmosphere in hydrostatic and ther- 
mal equilibrium confined between the radii r core and tl 
a roun d cores with mass larger than M crit , as explained in 

there are (at least) two 



3.1 



For masses below M cr a, 
solutions, corresponding to a low-mass and a high-mass 
envelope, respectively. 

In Figure [| we plot curves of total protoplanet mass M p i 
against core mass M core at different radii in a disk with 
a = 1CP 2 and M = 10~ 7 Mq yr -1 . In each frame, the 
different curves correspond to planetesimal accretion rates 
in the range 10 _11 -10 -6 M ffl yr -1 . The critical core mass 
is attained at the point where the curves start to loop 
backwards. 

When the core first begins to gravitationally bind some 
gas, the protoplanet is on the left on the lower branch 
of these curves. Assuming M core to be constant, as the 
core and the atmosphere grow in mass, the protoplanet 
moves along the lower branch up to the right, until the 
core reaches M cr i t . At that point the hydrostatic and 
thermal equilibrium approximation can no longer be used 
for the atmosphere, which begins to undergo very rapid 
contraction. Figure f| indicates that when the core mass 
reaches M cr u, the mass of the atmosphere is comparable 
to that of the core, in agreement with Bodcnheimer & Pol- 
lack (1986). Since the atmosphere in complete equilibrium 
is supported by the energy released by the planetesimals 
accreted onto the protoplanet, we expect the critical core 
mass to decrease as M core is reduced. This is indeed what 
we observe in Figure [|. 

For a = 1(T 2 and M = 1(T 7 M yr" 1 , the critical 
core mass at 5 AU varies between 16.2 and 1 M® as the 
planetesimal accretion rate varies between the largest and 
smallest value. The former result is in good agreement 
with that of Bodenheimer & Pollack (1986). Note that 
there is a tendency for the critical core masses to increase 
as the radial location moves inwards, the effect being most 
marked at small radii. At 1 AU, the critical mass varies 
from 17.5 to 3 M ffi as the accretion rate varies between the 
largest and smallest value, while at 0.05 AU these values 
increase still further to 42 and 9 M® , respectively. We note 
there has been some debate over the amount of grain opac- 
ity which should be used in these calculations. However, 
the results of Bodcnheimer & Pollack (1986) indicate that 
M cr i t does not depend sensitively on the grain opacity in 
the envelope. 

In Figure |5| we plot the critical core mass M cr n ver- 
sus the location r for three different steady disk models. 
These models have a = 10~ 2 and M = 10~ 7 M Q yr" 1 , 
= 10~ 2 and M = 10" 8 M Q yr" 1 , and a = 10~ 3 and 
yr -1 , respectively. 



a 

M = 



10~ s M P 



Here again, in each 



frame, the different curves correspond to planetesimal ac- 
cretion rates in the range 10 _11 -10 -6 M® yr -1 . Similar 
qualitative behavior is found for the three disk models, but 
the critical core masses are smaller for the models with 
M = M G yr" 1 , being reduced to 27 and 6 M® at 

0.05 AU for the highest and lowest accretion rate, respec- 
tively. 

These results indicate a relatively weak dependence on 
disk conditions except when rather high midplane tem- 
peratures T m > 1 , 000 K are attained, as in the inner re- 
gions. We indeed find similar values of M cr ;t for the three 
different models at radii larger than about 0.15-0.5 AU, 
where T m is lower than 1,000 K. Also, the fact that M cr n 
is similar for the two models with M — 10~ 8 M yr -1 
indicates that M cr n is more sensitive to the midplane 
temperature than to the midplane pressure. These two 
models have indeed similar T m , whereas P m varies signif- 
icantly from one model to the other. In the model with 
M = 10~ 7 M yr -1 , T m is significantly larger, hence the 
larger M cr n at small radii in this case. These results are 
consistent with the fact that M cr u depends on the bound- 
ary conditions only when a significant part of the envelope 
is convective (Wuchterl 1993), being larger for larger con- 
vective envelopes (Perri & Cameron 1974). In the model 
with a = 10~ 2 and M = 10~ 7 M yr -1 , we indeed find 
that the inner 70% in radius of the envelope is convec- 
tive at 0.05 AU, this value being reduced to 10% at 5 AU. 
When the envelope is mainly radiative, it converges rapidly 
to the radiative zero solution independently of its outer 
boundary conditions, so that M cr n hardly depends on the 
background temperature and pressure (Mizuno 1980). 

We note that it is unlikely there are planetesimals at the 
smallest radii considered here. Therefore, although criti- 
cal core masses for the same planetesimal accretion rates 
may be higher there, a lack of planetesimals may result in 
a fall in the core luminosity, making the critical core mass 
relatively small at these radii. 

4. PROTOPLANET MIGRATION AND PLANETESIMAL 
ACCRETION 

According to current models of planet formation 
(Safronov 1969; Wetherill & Stewart 1989), planetesimals 
form as the result of the coagulation of dust grains. Fur- 
ther accumulation through binary collisions then results in 
the formation of a core of several earth masses. After the 
critica l ma ss is attained, runaway gas accretion may begin 
(see § |3.ip . The core may form on a timescale of about 
10 5 yr at a distance of about 5 AU as a result of runaway 
planetesimals accretion (Lissauer & Stewart 1993). 

In addition, dynamical tidal interaction of a core of sev- 
eral earth masses with the surrounding disk matter be- 
comes important, leading to phenomena such as inward 
orbital migration and gap formation (Lin & Papaloizou 
1979a; Goldreich & Tremaine 1980; Lin & Papaloizou 
1993; Korycansky & Pollack 1993; Ward 1997b). For con- 
ditions under which the tidal interaction with the disk is 
linear, Ward (1986, 1997b) estimates inward migration (re- 
ferred to as type I) timescales t mig = —2r(dr/dt)~ 1 of 
about 2 x 10 5 (M p/ /M®) yr at 5 AU in a protoplane- 
tary disk similar to the minimum mass solar nebula. The 
timescales at other radii are expected to roughly scale as 
for the model disks considered here. They are also 
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somewhat shorter during the early phases of disk evolu- 
tion when the disk is more massive. 

However, for core masses of the magnitude we consider, 
the interaction may become nonlinear, leading to a reduc- 
tion in the inward migration (then referred to as type II) 
rate. To investigate the conditions for nonlinearity, Ko- 
rycansky & Papaloizou (1996) considered the perturbed 
disk flow around an embedded protoplanet assuming the 
disk viscosity to be negligible. They used a shearing-sheet 
approximation in which a patch, centered on the planet 
and corotating with its orbit, is considered in a 2D approx- 
imation. They found that the condition for nonlinearity 
or the formation of significant trailing shock waves in the 

1/3 

response is that r t (Q/c s ) > 0.5, where r t = r (M p i/M*) ' 
is a multiple of the Roche lobe radius. This condition ef- 
fectively compares the strength of the protoplanet 's grav- 
ity to local pressure forces. For the disk models consid- 
ered here, this condition is met for M p i = 10 M ffi for all 
disk mass accretion rates M at the inner radii. Even for 
M p i = 1 Mq, it is met at 0.1 AU for the higher accretion 
rates. Thus if we wish to consider core mass migration, 
nonlinear effects must be considered. These are expected 
to lead to a feedback reaction from the disk which cou- 
ples the orbital migration to the viscous evolution of the 
disk (Lin & Papaloizou 1986). However, this timescale 
can also be quite short, particularly for the models with 
higher accretion rates. For example, when a = 10 -3 and 
M = 10" 7 M vr" 1 , the viscous inflow timescale is on the 
order of 10 4 yr at 5 AU. 

The characteristic timescale of 10 4 — 10 5 yr obtained for 
migration and estimated for core formation suggests that 
cores of several earth masses form at about 5 AU and mi- 
grate inwards to small radii. In doing so, they continue to 
grow. As long as significant planetesimal accretion onto 
the core is maintained during the migration, the critical 
core mass M cr u remains significant and may not be at- 
tained before the core reaches small radii. At small radii, 
typically smaller than 0.1 AU, the planetesimal accretion 
rate decreases because the high temperatures have pre- 
vented the formation of planetcsimals there. The critical 
core mass is then also reduced below the actual core mass, 
so that runaway gas accretion can begin. Note though 
that the gas accretion phase may become longer tha n th e 



disk lifetime if the mass of the core is too small (see § 5T ) 
However, the build-up of a core massive enough through 
the merger of additional incoming cores may enable giant 
planets to form within the disk lifetime provided there is 
enough gas to supply the atmosphere. 

The accretion rate onto a protoplanet migrating in 
an approximately circular orbit through a planetesimal 
swarm with surface density S p can be estimated as: 



dM p 
dt 



2Y, p v R af. 



(22) 



Here we use a simple two dimensional model appropriate 
for a thin planetesimal disk from which accretion occurs 
onto a large protoplanet. The impact target radius of the 
protoplanet is a and vr is the relative velocity in colli- 
sions. We consider the case of near circular orbits. Then 
the relative velocity associated with a collision is expected 
to be equivalent to the Keplerian shear across the Roche 
lobe radius. Thus we adopt vr = 2rjSl as characteristic 



induced relative velocity. The factor / takes account of 
other effects such as gravitational focusing, which tends 
to increase the collision rate, and any local reduction in 
E p , which would be expected to occur if planetcsimals arc 
depleted locally and a gap tends to form (Tanaka & Ida 
1997). This might be expected for very slow migration 
rates but the total amount of accretion would be expected 
to be large in that case. 

We take a = O.Olrt as characteristic effective size of the 
protoplanet core, this representing the actual physical size 
of a core with density 3.2 g cm -3 at 0.75 AU. At smaller 
radii this is larger while at larger radii it is smaller. 

In this context we note that there is some uncertainty in 
the size of the target radius to be used because there may 
be a disk of bound planetesimals. Further, numerical tests 
indicate that the results we present below are not very sen- 
sitive to the magnitude of the target radius used because of 
the effects of gravitational focusing. This is also supported 
by the results of Kary, Lissauer & Grenzweig (1993). They 
considered the accretion of small planetesimals migrating 
inwards under the influence of gas drag by a protoplanet 
in fixed circular orbit (see below). Their results indicate 
that, for weak gas drag, which is appropriate for no reso- 
nant trapping, the impact probability typically differs from 
that assuming a target radius a = O.OItl by no more than 
a factor of about three as a varies between 0.0001r£ and 
tl . We comment that if the target radius was 6 times the 
actual core size, all accretion rates derived here would be 
underestimates for migration occurring with r < 5 AU. 

Thus for a simple estimate we use: 



dM p 
dt 



= 0.04/E p r£n. 



(23) 



We can estimate the total fraction of the total planetesimal 
mass accreted in a migration time t rniq to be: 



pi 



0.04/ 



7rS p r 2 dt 



U4J /r\L\- 
ir V r J 



nt r 



0.01 



Mp, 



2/3 
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mtg • 



(24) 



Characteristically, we find this fraction to be about 0.1 for 
/ = 1, M p i ~ 10 M ffi and flt mig — 10 4 , which are the 
expected characteristic values according to Ward (1997b). 
The expected fraction scales only weakly with protoplanet 

— 1/3 

mass, being proportional to M pl ' . But note too that 
the above may underestimate the accretion rate because 
larger relative velocities may be induced if there are mul- 
tiple close scatterings. 

It is of interest to compare the accretion rate expected 
from the above two dimensional model with predictions 
based on the standard accretion formula with gravitational 
focusing for three dimensions (Lissauer & Stewart 1993). 
This gives: 



dMpi 
dt 



oVR 



2h r 



2GM p i \ 

aV R J 



(25) 



where h p = vr/Q, is the semithickness of the planetesimal 
distribution. Using the same estimate for vr as above, this 
gives (assuming the dominance of the second gravitational 
focusing term in the brackets): 
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(26) 



For the same parameters as used above this gives the same 
prediction as equation (^) with / = 2. Both these expres- 
sions and our simulations give consistent results suggesting 
that the migrating protoplanet accretes as if it were in a 
homogeneous medium without a gap forming in the plan- 
etesimal distribution. Whether such a gap forms should be 
reliably determined by the two dimensional calculations. 

Given that the disk is expected to contain at lea st a bout 
8 M ffi within 5 AU in the early stages (see § 2J3 and 
Table |]), these estimates indicate that an accretion rate 
Mcore of at least about 1CP 6 M ffl yr -1 is likely to be main- 
tained during orbital migration in the present disk models. 
Efficient gas accretion is then unlikely to start until small 
radii are reached, at least in the early phases of the disk 
lifetime. Note too that the fractional accretion rate given 
by equation (p4j) is not expected to increase indefinitely 
with t m i g because of the tendency to form a gap (Tanaka 
& Ida 1997) which would then be expected to cause a re- 
duction in /. 

In order to verify the above conclusions, we have per- 
formed simulations of migrating protoplanets with M v i 

1 M ffi and migration times, t m i g , measured at 1 AU, of 
between 2 x 10 3 yr and 10 4 yr. We have also considered the 
case M p i = 1 M ffi and t m i g between 2 x 10 4 yr and 10 5 yr. 
The protoplanet was taken to be in a quasi-circular or- 
bit with ln(r) decreasing on the specified timescale. The 
protoplanet was assumed to start at 1 AU but the results 
can be scaled to any other initial radius in the usual way. 
The 256 planetesimals were initially regularly spaced be- 
tween 0.6 AU and 0.8 AU, as indicated in Figure ^. The 
protoplanet was allowed to migrate through them. To esti- 
mate the fraction of planetesimals accreted, just as above, 
we assumed that any particle approaching the protoplanet 
within 0.01 tl was accreted. 

We remark that the migration was imposed here, be- 
ing possibly due to interactions with the disk. However, 
a migration mechanism based on the scattering of plan- 
etesimals alone has been noted by Murray et al. (1998). 
Their mechanism requires the protoplanet orbit to have 
a non zero eccentricity, but this could be damped by the 
disk interaction (Artymowicz 1994). 

For the protoplanet masses and migration rates we con- 
sider, we found significant accretion of planetesimals. We 
here present two examples. The final distribution of the 
planetesimals after a protoplanet of 1 Mq has migrated 
through them with t m i g = 2 x 10 yr is indicated in Fig- 
ure ^. In this case, about 24% of the planetesimals initially 
present were accreted. Even more planetesimals were ac- 
creted with slower migration rates. 

The final distribution of the planetesimals after a proto- 
planet of 10 M ffi has migrated through them with t m i g = 

2 x 10 3 yr is given in Figure ||. In this case, about 23% 
of the planetesimals initially present were accreted. Given 
that the disk models typically contain at least about 8 M ffi 
interior to 5 AU, we conclude that enough accretion occurs 
during the migration to prevent gas accretion as long as 
planetesimals are present. 

It is of interest to compare the results of our simulations 
with those of Kary et al. (1993). These authors consid- 



ered the accretion of small bodies migrating inwards under 
the influence of gas drag onto a protoplanet in fixed cir- 
cular orbit. Although this is not an identical situation to 
the one we consider here, it is similar enough to make a 
comparison interesting. 

An important aspect of these simulations is that gas 
drag causes eccentricity damping as well as inward mi- 
gration. The eccentricity damping allows particles to 
be trapped in resonances such that they stop migrating, 
maintaining near circular orbits, with the consequence 
that close approaches to the protoplanet are avoided. The 
resonant interaction causes eccentricity growth at a rate 
governed by the migration rate. This is because energy 
and angular momentum transfer from the protoplanet is 
in the wrong ratio for keeping the other body in circular 
orbit. Without eccentricity damping, near circular orbits 
for the particles could not be maintained (see, for example, 
Lin & Papaloizou 1979b for a discussion). 

To relate the migration and damping rates, we consider 
a particle with semi-major axis a p and eccentricity e un- 
dergoing resonant interaction with a protoplanet in a near 
circular orbit with semi-major axis a p i. Under these con- 
ditions, the Jacobi integral: 



Jc 



Op 

dpi 



3/2 



(27) 



is conserved for the particle. Thus, changes to a p and e 
are related by: 



dt 



3/2 



y/l - e 2 



v/I^-^. (28) 
a p at 



When an inwardly migrating protoplanet pushes an inte- 
rior particle in front of it maintaining a fixed period ratio, 
just as for the protoplanet, — 2a p (da p /dt)~ 1 — t rn i g . Then 
equation d28| ) implies that the eccentricity increases. If the 
orbit is to maintain a finite eccentricity, this rate of in- 
crease has to be balanced by the the eccentricity damping 
or circularization rate due to dissipative processes such as 
gas drag. If the circularization time is t c i rc — —e^de/dt)^ 1 , 
an equilibrium eccentricity can be maintained, such that 
for small e: 



3/2 

"7" \ _ 1 



a 



tr 



(29) 



Equation ([29|) also applies to the case of a fixed proto- 
planet orbit and particles migrating inwards due to gas 
drag. In that case a p i < a p , and the resonant interaction 
gives positive rates of increase for both a p and e. These 
are balanced by the inward migration rate due to gas drag 
and the orbital circularization rate respectively. 

The process of resonant trapping is accordingly expected 
to be similar in the cases of a free particle with inward 
migrating protoplanet and particle migrating inwards to- 
wards a protoplanet on fixed circular orbit. But it is im- 
portant to note that the physical processes causing migra- 
tion and circularization may both be very different. 

An interpolation of the results of Kary et al. (1993) in- 
dicates that resonant trapping is important for t m i g > 6 x 
10 4 yr for a 1 M ffi protoplanet at 1 AU, and t mig > 10 4 yr 
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for a 10 M ffi protoplanet at 1 AU. These migration times 
are longer than those proposed by Ward (1997b), or the 
disk viscous timescale at 1-5 AU in the early stages of the 
protoplanetary disk considered here. The expectation is 
that resonant trapping will not be important then. 

When there is no resonant trapping, the fraction of bod- 
ies accreted is similar to that found here, ranging between 
10 and 40 percent. We comment that the high rate of 
eccentricity damping required for effective resonant trap- 
ping is only likely to be obtained for small bodies. From 
Kary et al. (1993), when marginal trapping occurs for 
a 10 M ffi protoplanet approached by small bodies with 
tnug = 10 4 yr, t circ ~ 30 yr at 1 AU. As t circ is propor- 
tional to the radius of the body, values less than 30 yr 
require the radius of the body to be smaller than 20 m at 
1 AU. But note that the relative effectiveness of gas drag 
tends to increase at smaller radii. Thus resonant trapping, 
should it occur, is more likely in the inner regions of the 
disk. 

Equation (|?|) suggests that t m i g and t circ should scale 
together for marginal resonant trapping occurring with the 
same orbital configuration. Thus a 10 M ffi protoplanet at 
1 AU migrating inwards with t m i g — 2 x 10 5 yr should res- 
onantly trap bodies, causing them also to migrate inwards 
maintaining a fixed period ratio, if t C irc < 6 x 10 2 yr. We 
have verified that trapping occurs when t C i rc — 10 2 yr. 

However, it appears that resonant trapping as a result 
of gas drag is unlikely for planetcsimals with radii larger 
than about 10 km and the accretion rates should then be 
similar to those found here. Departures are to be expected 
only for small bodies at the slowest migration rates. 

5. DISCUSSION AND SUMMARY 

5.1. Formation of Short-Period Giant Planets 

The above suggests that a protoplanetary core formed at 
about 5 AU which migrates inwards will not attain the crit- 
ical core mass, above which runaway gas accretion starts, 
before it reaches small radii ~ 0.05-0.1 AU where plan- 
etesimals no longer exist. Runaway gas accretion onto a 
small core can then occur at these radii. However, if the 
core is too small, the gas accretion phase may be longer 
than the disk lifetime. 

Even for core masses in the range 15-20 M ffi , the build-up 
of a massive atmosphere may take a time ~ 10 6 yr (Boden- 
hcimer et al. 1998). The reason for this is that once the 
core starts to accrete a significant atmosphere, energy pro- 
duction occurs through its gravitational contraction. The 
luminosity produced then slows down the evolution. An 
estimate of the evolutionary timescale at this stage can be 
obtained by noting that the luminosity should be equiva- 
lent to that required to make the core critical (and hence 
to enable runaway gas accretion to begin), assuming it to 
be produced by planetesimal accretion. We may thus esti- 
mate the time scale as the Kelvin-Helmholtz time for our 
already calculated critical core mass models. This is given 
by: 

tKH = \E\/L core , (30) 

where E is the total internal and gravitational energy of 
the gas. 

We have calculated txH in this way for critical cores at 
different radial locations in a disk model with a — 10 -2 



and M = 10 7 M Q yr 1 . The results are presented in 
Table |. Typically, we find that t KH ~ 10 6 - 10 7 yr for 
core masses in the range 10-20 M ffi for radii larger than 
0.075 AU. The core masses required to get such a char- 
acteristic timescale increase rapidly interior to 0.06 AU. 
However, we note that they decrease as the mass transfer 
rate through the disk does. The fact that fairly large core 
masses are required to give evolutionary timescales com- 
parable or less than the expected disk lifetime means that 
mergers of additional incoming cores may be required in 
order to produce a core of sufficient mass that real runaway 
gas accretion may begin. 

Table [l] indicates that, in a disk with a = 10~ 3 and 
M = 10~ 6 or 10" 7 M Q yr" 1 , 40 M e of planetesimals are 
contained within 1 or 5 AU. Therefore, the timescale for 
building-up a core with a mass between 20 and 40 M ffi 
at small radii is typically the timescale it takes for plan- 
etesimals to migrate from 1 or 5 AU down to these small 
radii. According to Ward (1997b), the migration timescale 
of cores of a few tenths of an earth mass located at 1 or 
5 AU is at most 10 6 yr in such a disk, and it decreases 
with increasing core mass. Therefore, if planetesimals can 
be assembled into cores of at least a few tenths of an earth 
mass at these radii on a reasonably short timescale, a mas- 
sive core could be obtained at small radii on a timescale 
much shorter than for in situ formation. 
If the disk has a = 10~ 2 , 40 of planetesimals are con- 
tained within 5 or 11 AU. In this case again, the migration 
timescale of cores of a few tenths of an earth mass located 
at 5 or 11 AU is about 10 6 yr, so that the above discussion 
still holds. 

A massive core can be built-up through the merger of 
additional incoming cores either after having stopped at 
small radii or on its way down to small radii (where it 
would still be expected to be stopped). The former pro- 
cess resembles that discussed by Ward (1997a). The latter 
scenario would occur if more massive cores, which migrate 
faster, overtake less massive cores on their way down. 

Supposing that a protoplanetary core massive enough 
can be built-up on its way down to small radii and that it 
continues to rapidly move inward until it gets interior to 
the disk inner boundary, it can only accrete the gas which 
is in its vicinity, i.e. typically the amount of gas con- 
tained within ~ 0.05-0.1 AU. Since the core is expected 
to reach these radii early in the life of the protoplanetary 
disk, ther e may still be an adequate amount of gas there 
(see § ||~J|) for it to build-up a large envelope and become 
a giant planet. However the conditions for that to happen 
are rather marginal. 

If the protoplanetary core is stopped at some small ra- 
dius before the disk is terminated, it might be able to re- 
tain contact with disk gas. In that case it might be able to 
accrete enough gas supplied from the outer disk by viscous 
evolution to build-up a massive atmosphere. 

The question arises as to the nature of any process that 
can halt the migration. One might suppose that because 
the viscous evolution timescale of the disk increases as 
the disk gets older, protoplanet migration, if type II, gets 
slower and slower, and may halt altogether when the disk 
dissipates (e.g., Trilling et al. 1998). However, given that 
the migration timescale is shorter at smaller radii, if no 
mechanism halts it there, very fine 'tuning' would be re- 
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quired to produce the large fraction of extrasolar planets 
found on very close orbits in this way. 

The disk may be terminated by a magnetospheric cavity 
(Lin et. al 1996) such that tidal torques producing inward 
migration vanish once the protoplanet enters it. However, 
as a result of such an entry, contact with the disk is lost 
and further accretion of gas may be difficult. Also magne- 
tospheric cavities may not extend up to a few tenths of an 
AU, where some of the extrasolar planets have been found 
to orbit. 

The way in which the migration of a protoplanet would 
be halted in the inner regions is not yet clear (e.g., Bo- 
dcnhcimer et al. 1998). In this context, we remark that 
migration rates are usually considered in relation to a stan- 
dard a disk model in which the disk midplane tempera- 
ture increases inwards. In this situation inward migration 
is expected in general (Ward 1986). However, if the in- 
ner disk is terminated through interaction with a stellar 
magnetic field, physical conditions may start to differ in 
the interaction zone where magnetic field lines penetrate 
the disk. Additional energy and angular momentum trans- 
port mechanisms due to a wind for exa mpl e may start to 
become important (see discussion in § 2.3). As a result, 
an inward midplane temperature decrease might be pro- 
duced. It may then be possible that migration could be 
halted such that the protoplanet retains contact with disk 
gas. 

In the context of several cores interacting together, we 
note that stellar tides are unlikely to provide enough ec- 
centricity damping for resonant trapping to occur for the 
migration rates we consider. The circularization time re- 
sulting from stellar tides is given by (Goldreich and Soter 
1966): 



icirc (yr) ~ 2.8 x 1Q- 5 Q 



M p i P ( a p i 
M7l day \~Rp~ 



Here P Q denotes the orbital period, R p the radius of the 
planet and Q is the usual tidal Q-value. For a 10 earth 
masses planet orbiting a solar mass at 0.05 AU, one typi- 
cally finds i c i rc ~ 3 x 10 4 Q yr. As Q > 1, this must exceed 
the migration times considered here, and so resonant trap- 
ping is unlikely. This is because normally t clrc -C t m i g is 
required (Kary et al. 1993, Lin & Papaloizou 1979b and 
see eq. p9|). Similar conclusions can be obtained if stellar 
tides acting on a Jovian planet are considered (Lin et al. 
1998). 

However, conditions may be different later in the life 
of the disk when the viscous evolution time is longer and 
type II migration rates are slower. A model of the proto- 
stellar disk in which giant planets form at about 5 AU at a 
later stage in the life of the disk after about 10 6 — 10 7 yr, 
when the viscous time is of comparable magnitude, has 
been considered by Trilling et al. (1998). In this situa- 
tion the planet is able to open a gap and undergo inward 
type II migration on the viscous timescale. These authors 
consider outward torques due to tidal interaction with the 
central star and Roche lobe overflow as well as torques due 
to disk interaction. Assuming disk dispersal on a similar 
timescale they are able to produce giant planets at a range 
of orbital radii. These may then undergo orbital instabil- 
ity, leaving one inner planet and one or several partners at 
larger radii (e.g., Weidenschilling & Marzari 1996; Rasio 



& Ford 1996). In this regard, it is of interest to note that 
the only system detected so far which may be multiple, 
55 Cnc (Butler et al. 1997), has a planet at 0.11 AU and 
maybe another planet beyond 4 AU. 

The above scenario thus might be able to produce short 
period planets in the late stages of the life of the disk. In 
contrast, the processes which are the focus of this paper 
result in the short period planets originating early in the 
life of the disk. They would more likely result in a single 
planet at ~ 0.05 — 0.1 AU than at intermediate radii and 
do not necessarily produce other giant planets at larger 
radii as a result of gravitational scattering processes. In 
this respect, the outcome would be in good agreement with 
the observations to date. 

5.2. Summary 

In this paper we have investigated how the critical core 
mass associated with a solid protoplanet varies with loca- 
tion in a protoplanetary disk. In the past, work has con- 
centrated on in situ formation, mainly at orbital distances 
corresponding to the neighborhood of Jupiter. However, 
the importance of orbital migration has recently been indi- 
cated by both observations (see Marcy & Butler 1998 and 
references therein) and theory (see Lin & Papaloizou 1993 
and references therein; Ward 1997b). This suggests that 
the behavior of the critical core mass as a function of loca- 
tion in and mass transfer rate through the protoplanetary 
disk should be considered. 

We constructed steady state protostellar disk models 
with constant values of a = 10~ 2 and a = 10 -3 . A 
range of accretion rates varying between M = 10 -6 and 
10~ 9 M yr -1 were considered. We calculated analytic 
piece-wise power law fits to the curves (v) (£) obtained 
from numerical calculations. These fits can be used to 
solve the diffusion equation governing the disk evolution 
for a wide range of disk parameters. 

We constructed protoplanet models with a solid core 
and a gaseous envelope in hydrostatic and thermal equi- 
librium. We calculated the critical core mass, M cr u, above 
which the atmosphere cannot remain in complete equilib- 
rium but must begin to undergo a very rapid contraction. 
Where they can be compared, these critical core masses 
agree with the crossover masses obtained from more so- 
phisticated time dependent calculations (Bodenheimer & 
Pollack 1986, Pollack et al. 1996). 

We found that, for a fixed core accretion rate, M cr u typ- 
ically increases by factors of 2-3 between 5 and 0.05 AU for 
disk parameters believed to be typical of the early stages 
of the disk evolution (i.e. a — 10 -2 and a gas accretion 
rate M = 10~ 7 Mq yr -1 ). For a core accretion rate of 

M core = 10 -6 M ffi yr -1 , the critical core mass is found 
to be about 16 M ffi at 5 AU, in agreement with Boden- 
heimer & Pollack (1986), and it decreases with M core . At 
radii smaller than about 0.05-0.1 AU, local conditions may 
not enable the planetesimals, required to produce the ac- 
cretion luminosity to support the atmosphere, to exist. 
Therefore the critical core mass is reduced at these radii. 

We considered these results in the context of the mi- 
gration of protoplanetary cores with mass in the range 
1-10 M$. According to recent estimates of the effects 
of tidal interactions with the disk (Ward 1997b), the mi- 
gration timescale for such cores may be between 10 and 
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10 5 yr at 5 AU, being comparable to their proposed forma- 
tion time. In order to investigate whether the planetesimal 
accretion would continue under migration, we performed 
simulations which indicate that a protoplanet migrating 
at the proposed rate is likely to accrete 24% or more of 
the planctesimals initially interior to its orbit. Thus ac- 
cretion is likely to maintain the core luminosity such that 
attainment of the critical mass does not occur until small 
radii ~ 0.1 AU are reached, where planetesimals no longer 
exist. 

Although runaway gas accretion can then begin onto 
small mass cores at these small radii, the timescale for 
building-up a massive envelope becomes longer than the 
disk lifetime if the core is too small. However, cores mas- 
sive enough can be built-up through mergers of additional 
incoming cores on a timescale shorter than for in situ for- 



mation. 

The above considerations can lead to the preferential 
formation of short-period planets with semi-major axis 
less than about 0.1 AU, on a timescale shorter than that 
derived for in situ planet formation. 
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APPENDIX 

A. BOUNDARY CONDITIONS FOR THE DISK VERTICAL STRUCTURE 



We derive here the surface pressure P s and temperature T s used to compute the disk vertical structure (see § [2.1.2 ). 
To get P s we rewrite equation ([!]) under the form: 

dp n 2 z , k s 

~r = , (ai) 

where r(z) = L pndz is the optical depth. We then integrate this equation over r between the surface of the disk and 
infinity, where the pressure is zero. This leads to: 

P s = / —dr. (A2) 

Jt(H) k 

We define the disk surface such that the atmosphere above the disk is isothermal. The mass density above the disk then 
varies like exp[— Q 2 {z 2 — H 2 ) I (2cf)l , where c s is the (constant) sound speed in the atmosphere. The main contribution 



to the integral in equation (A2) comes from values of z starting from H and extending over a range of a few c s /f2 which is 
significantly less than H. Thus the integral can be evaluated with sufficient accuracy by taking fl 2 z and k to be constant 
and equal to their values at the disk surface. This leads to: 

P, - (A3, 

where T a b = J H dr is the optical depth above the disk. 

We now calculate the surface temperature T s . The radiative flux at the disk surface can be written under the form: 

F S = F+-F_ 1 (A4) 

where F + and F_ are the radiative vertical fluxes at the surface of the disk directed respectively along the positive and 
negative z. In other words, F + and Fl come respectively from inside and above the disk. The flux F_ can be w ritte n in 



term of the background temperature T}, in which the disk is embedded as F_ = uT*. Using equations and (A4), we 
then get: 

F+ = ^-M st il 2 + oT£. (A5) 
An other expression of F + can be obtained by using the fact that the energy density at the disk surface is given by: 

E S = -(F++F_), (A6) 
c 

together with the energy equation: 

V • F = pne (aT 4 - E) , (A7) 

where F is the vector representing the flux of radiative energy, E is the energy density, c is the speed of light and a = Aa/c 
is the radiation constant. In the thin disk approximation, the temperature gradient in the vertical direction is much larger 
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than that in the horizontal dire ction , so t hat V • F ~ dF/dz. Using the expression (Q) and the relation between F_ and 
Tj,, we then get from equations (|A6|) and (A7) written at the disk surface: 



2aT? - oTi 



9a (c 2 s ) s n 



8k s 



(A8) 



where v has been expressed in term of a. By comparing equations ( A5) and ( A8) we obtain the following equation for T s : 

(A9) 



9a (ci) ft 3 . 
2a (Tt - Th) ^ ^-Msttf = 0. 



8k., 



8tt 



B. ANALYTICAL FITS OF THE VERTICAL STRUCTURE MODELS 

To compute the evolution of a non-steady a-disk, it is necessary to solve the diffusion equation (||). To do this, 
M st = 37r(^)E has to be specified as a function of E at each radius. We found it convenient to have an analytic fit to the 
curves M st (E) . 

We utilize fits in which these curves are approximated by three different power laws corresponding respectively to the 
optically thin, intermediate and optically thick regimes. The index of these power laws is independent of the radius r and 
the parameter a. However, the multiplicative constant characterizing each of them does vary with both r and a. We give 
this dependence below. 

Figure ^ shows a schematic plot of the fits M st (E). We have represented log 10 (^M s tj vs. log 10 (S) with an arbitrary 

scale at two different arbitrary radii n and r 2 such that r 2 > r%. When the radius is increased, the point at the 
transition between the optically thin and intermediate regimes on this logarithmic representation moves up along the 
straight line with equation is y = 3.1a; + c'. We now give details of these power laws, which give a good fit (see below) 
for M s t < 1CU 4 MQ/year and 1CU 5 < a < 10" 1 . In the following expressions, M st , r and E are in cgs units and the 
logarithms are to base 10. 

In the optically thin regime, as long as E < Ei with log(Ei) = (ci — c')/2.1, or equivalently M st < M s t,i with 
log (M st ,i) = (3.1ci - c')/2.1, we have: 

log (M st ) = ct + log (E) . (Bl) 

In the intermediate regime, for Ei < E < E2 with log(E2) = (03 — C2) /0.9, or equivalently M stj i < M st < M s t,2 with 
log (Af st , 2 ) = (2c 3 - l.lc 2 )/0.9, the fit is: 

log(M st ) =c 2 + 21og(E). (B2) 
Finally in the optically thick regime, for E > E 2 , or equivalently M st > Af st 2 , we use: 

log(M flt ) =c 3 + l.llog(E). (B3) 

The parameters ciand C3 are related to r and a in the following way: 

log(ci) = 0.9360636 + 0.1195816 log(a) + [0.0233002 - 0.0061733 log(a)] log(r), (B4) 

log(c 3 ) = 0.7782080 + 0.0545617 log(a) + [0.0366565 - 0.0019087 log(a)] log(r), (B5) 
whereas c' depends only on a: 

c = 16.0897161 + 2.0665 log(a), (B6) 

and 

C2 = —2J-- {B7) 

We note that these fits can be used either to compute E from M s t or M s t from E. 

We calculate the error corresponding to these fits by solving the vertical structure (i.e. calculating E) at 50 different 

radii from 0.01 to 100 AU and for 50 values of M st between 10~ 10 and 10 Mq yr _1 (both the values of r and M st are 
equally logarithmically spaced). We then recalculate E using the fits and M st as an input parameter. We find that the 
average error is 24, 17, 14, 10 and 10%, respectively, for a = 10~ 5 , 10~ 4 , 10 -3 , 10 -2 and 10 _1 . The maximum error is 
between 43 and 50% for these values of a. If alternatively we recalculate M st using the fits and E as an input parameter, 
the average error is 36, 22, 18, 13 and 16% whereas the maximum error is 107, 55, 48, 42 and 103% for the same values 
of a from 10~ 5 to 10 _1 . To solve the radial diffusion equation, we need to calculate M st from E (see § JT^). We see that, 
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except for a = 10 5 and 10 , the fits give a good approximation. For a = 10 \ we can get the maximum and average 
errors down to about 50 and 10%, respectively, by limiting the calculations to M st > 10~ 8 Mq yr _1 . For a = 10~ 5 , the 
maximum error gets down to 76% as M st varies between 10 -8 and 3 x 10~ 5 Mq yr _1 , while the average error does not 
change. 
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Table 1 

Mass of planetesimals contained within a radius r. 
Listed are the radius r of the disk location in AU, a, the gas accretion rate M through the disk in Mq yr -1 
and the estimate for the mass of planetesimals contained within the radius r GIVEN BY M p (r) = 1CT 2 X 7rr 2 E(r) IN 

Me. 



r 

(AU) 


a 


M 

(M Q yr- 1 ) 


M p (r) = IO' 2 x 7rr 2 S(r) 
(M e ) 


1 


10~' 2 


io- b 


3.5 






io- 7 


0.6 


5 




1(T 6 


22.8 






io- 7 


7.6 


1 


1Q- 3 


10~ 6 


30.1 






io- 7 


4.0 


5 




IO" 6 


138.5 






io- 7 


41.0 



Table 2 

Kelvin-Helmholtz timescale for different critical core masses. 
The first column gives the radius r OF the disk location in AU, THE second column gives the internally 

GENERATED PROTOPLANET LUMINOSITY AS WHAT WOULD BE DERIVED FROM A CORE ACCRETION RATE Mcore IN M® YR" 1 , AND 
THE THIRD AND FOURTH COLUMNS GIVE THE CORE MASS M core AND TOTAL MASS M pl IN M® , RESPECTIVELY. THE FIFTH 
COLUMN GIVES THE KELVIN-HELMHOLTZ TIME t K H IN YR. THE CALCULATIONS WERE PERFORMED FOR A DISK MODEL WITH 

a = 1CT 2 and M = 1CT 7 Mn YR -1 . 



r 






M p i 


tKH 


(AU) 


(M e yr" 1 ) 


(Me) 


(Me) 


(yr) 


0.05 


IO"" 


42.0 


56.9 


4.7xl0 e 




io- 7 


28.0 


32.7 


3.0xl0 7 


0.075 


10~ 6 


21.5 


28.2 


1.9xl0 6 




io- 7 


15.0 


18.9 


1.4xl0 7 


0.10 


10~ 6 


19.5 


24.8 


1.4xl0 6 




io- 7 


14.0 
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Fig. 1. — M s t in units Mq yr~ 1 vs. E in units g cm~ 2 using a logarithmic scale for a = 10 -2 (a) and 10~ 3 (b). Both the curves 
corresponding to the numerical calculations (solid line) and the fits (dashed line) are shown. The label on the curves represents the radius, 
which varies between 0.01 and 100 AU. 
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Fig. 2. — Shown is H/r (a), S in units g cm 2 (b) and T m in K (c) vs. r in units AU using a logarithmic scale. In each plot, the different 
curves correspond to M s t = (solid line), 10~ 7 (dotted line), 10 — 8 (short-dashed line) and 10 — 9 (long-dashed line) Mq yr -1 . Here 

a = 10- 2 . 
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Fig. 3. — Same as Figured but for a = 10 3 
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Fig. 4. — Plots of total mass, M p i in units M9, vs. core mass, M core in units M9, at different locations r in a steady state disk model with 
a = 10~ 2 and gas accretion rate M = 10~ 7 Mq yr -1 . From left to right and top to bottom, the frames correspond to r = 0.05, 0.06, 0.15, 0.5, 
1 and 5 AU, respectively. The midplane temperature and pressure at these locations are indicated above each frame. Each frame contains six 
curves which, moving from left to right, correspond to core luminosities derived from planetesimal accretion rates of M core = 10 -11 , 10~ 1() , 
10~ 9 , 10~ 8 , 10~ 7 and 10~ 6 M ffi yr" 1 , respectively. The critical core mass is attained when the curves first begin to loop backwards when 
moving from left to right. 
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Fig. 5. — Critical core mass M cr n in units Mgj vs. location r in units AU in a steady state disk model with viscous parameter a and gas 
accretion rate M. The different plots correspond to a = 1CT 2 and M = 1CT 7 M yr~ x (a), a = 1CT 2 and M = 10~ 8 M Q yr _1 (b), and 
a = 10 -3 and M = 10~ 8 Mq yr _1 (c). Each plot contains six curves which, moving from bottom to top, correspond to core luminosities 
derived from planetesimal accretion rates of M core = 10~ n , 10~ 10 , 10~ 9 , 10~ 8 , 10~ 7 and 10 -6 Mq yr _1 , respectively. 
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Fig. 6. — This shows the initial distribution of planetesimals (represented by crosses) between 0.8 AU and 0.6 AU in the x-y plane in units 
AU. The initial circular orbit of the protoplanct at 1 AU is also indicated (solid line). The protoplanet is represented by the open square. 
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Fig. 7. — The final distribution of planctcsimals (represented by crosses), in the x-y plane in units AU, after a protoplanet of mass 1 Mj 
has migrated inwards with t m i g = 2 X 10 4 yr. The initial protoplanet orbit is indicated by the outer circle (dashed line) and the final one by 
the inner circle (solid line). The protoplanet is represented by the open square. About 24% of the planetesimals were accreted. 
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Fig. 8. — The final distribution of planetcsimals (represented by crosses), in the x—y plane in units AU, after a protoplanct of mass 10 
has migrated inwards with t m i„ = 2 X 10 3 yr. The initial protoplanet orbit is indicated by the outer circle (dashed line) and the final one by 
the very small inner circle (solid line). The protoplanet is represented by the open square. About 23% of the planetcsimals were accreted. 
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Fig. 9. — Schematic plot of the fits M st (S). The solid lines represent a fit of the curves log 10 (-M 3 t) vs. log 10 (£) with an arbitrary scale at 
two different arbitrary radii ri and T2 such that T2 > ri . The dotted line indicates the separation between the optically thin and intermediate 
regimes. 



